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ABSTRACT 

In structural analysis the amount of computational time necessary for a solu- 
tion is proportional to the number of degrees of freedom times the bandwidth 
squared. In implicit time analysis, this solution must be calculated at each 
discrete point in time. If, in addition, the problem is nonlinear, then this 
solution must be iterated at each point in time. If the bandwidth is large, 
the size of the problem that can be analyzed is severely limited. 

The multigrid method is a possible algorithm that can make this solution much 
more computationally efficient. This method has been used for years in compu- 
tational fluid mechanics. It works on the fact that relaxation is very effi- 
cient on the high-frequency components of the solution (nearest-neighbor 
interactions) but is not very efficient on the low-frequency components of the 
solution (far interactions). The multigrid method relaxes the solution on a 
particular model until the residual stops changing, which indicates that the 
solution contains the higher frequency components. A coarse model is then gen- 
erated and relaxed for the lower frequency components of the solution. These 
lower frequency components are then interpolated to the fine model. 

In computational fluid mechanics the equations are usually expressed as finite 
differences. A coarse model is generated by just doubling the grid size and 
using a Green's integral theorem to obtain the forcing function on the coarse 
grid. Linear interpolation is used to transfer the lower frequency solution 
back to the fine grid. 

In structural dynamics the equations are usually expressed as finite elements. 
Neighbor elements need not be connected. The process of condensing a fine 
model into a coarse model and interpolating the low-frequency solution to the 
fine model will be studied in this work. 


OBJECTIVE 


The objective of this work is to use an implicit time march solution to study 
nonlinear structural dynamics. The work will be done in three phases. The 
first phase, a beam structure, will have application in a multishaft, combined 
lateral, torsional, and axial rotordynamic analysis. The second phase, a plate 
structure, will have application in bladed disk vibration with coulomb damping. 
The third phase, a full three-dimensional structure, will have application in 
space structures. 

To aid the reader, a symbols list has been included in the appendix. 


IMPLICIT TIME MARCH SOLUTION OF NONLINEAR STRUCTURAL DYNAMICS 

• BEAM-MULTISHAFT, COMBINED LATERAL, TORSIONAL, AND AXIAL ANALYSIS 

• PLATE-BLADE VIBRATION WITH COULOMB DAMPING 

• THREE-DIMENSIONAL SPACE STRUCTURES ANALYSIS 
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NUMERICAL INTEGRATION 


The numerical integration method is based on a Nordsieck-like method. The 
displacement, velocity, and acceleration are defined at an initial time. A 
modified Taylor series is used to calculate the displacement, velocity, and 
acceleration at the advanced time. The Lagrange remainder term, the time 
derivative of the acceleration, is calculated from the equations of motion 
at the advanced time. The constants a and |J are determined so that the 
method is stable as time approaches infinity. 

This method of integration for a first-order differential equation is Gear’s 
method (Gear, 1971). Zeleznik showed that this method could be used on higher 
order equations (private communication with F.J. Zeleznik at NASA Lewis 
Research Center in 1979). Kascak (1980) showed that for a third-order inte- 
grator used on a linear second-order differential equation the method is 
unconditionally stable. 


LET R(t) BEAN n ELEMENT VECTOR OF NODAL DISPLACEMENTS AND 
V(t) * R A(t> = V 
MODIFIED TAYLOR SERIES 

R(t) = RIO) + V(0)t + 1 A(0)t 2 + 1 oA(E)t 3 
2 6 

Vlt) =V(0) + A(0)t + lpAIE)t 2 
Alt) * A(0) + A(E)t 


WHERE 


0<E <t 

AND a AND 0 ARE DETERMINED SO THAT THE METHOD IS NUMERI- 
CALLY STABLE AS t-» 
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NUMERICAL STABILITY 


The numerical stability of the integration method can be examined by substi- 
tuting the displacement, velocity, and acceleration into the linear equations 
of motion, and solving for the time derivative of the acceleration. As time 
approaches infinity the dominate term on each side of the equation has the 
stiffness matrix as a premultiplier. The time derivative of the acceleration 
is proportional to the initial acceleration divided by the time. If the time 
derivative is substituted into the modified Taylor series and if a is set to 
3 and |3 is set to 2, then the acceleration is zero and the velocity is con- 
stant. The eigenvalues become zero and one. 


MA + CV + KR = F 

(tM + 1 pt 2 C + I ot 3 K ) A(U - F - + tC + I t 2 l<) A(0) - (C + tK)V(O) - KR(0) 

AS t- °° 

R - R(0) + V(0)t 
V = V(0) + ^1 - |(1)) AW 

A = (1- (lV\ A(0) 

\ \ a// 

LET a * 3 AND p ■ 2 

R * R(0) + V(0)t V*V(0) A = 0 
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ITERATIVE SOLUTION 


If the initial displacement, velocity, and acceleration, and an initial esti- 
mate of the time derivative of the acceleration are given, then estimates of 
the advanced displacement, velocity, and acceleration are given by using the 
modified Taylor series. If a correction to the estimate of the time derivative 
of the acceleration is given, then new estimates of the displacement, velocity, 
and acceleration are given by the modified Taylor series. The correction to 
the time derivative of the acceleration can be found from the equations of 
motion. 


GIVEN R(0), V(0), A(0), AND A(E;)~A (0) 

THEN R (0) = R(0) + V(0)t + I A(0)t 2 + I aA (0) t 3 

2 2 

V<°> = V(0) + A(0)t + 1 pA (0) t 2 

A (0) = A(0) + A (0) t 

LET A(^)=A (0) + AA 

THEN R(t) = R (0) + 1 aAAt 3 

6 

V(t) = V (0> + 1 pAAt 2 
Alt) = A (0) + AAt 
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NONLINEAR EQUATIONS OF MOTION 


The nonlinear equations of motion are the sum of both the static and dynamic 
forces for each element. As such, the equations are functions of the displace- 
ment, velocity, acceleration, and time. If the modified Taylor series is 
substituted into the equations of motion using the iterative form, then the 
equations of motion become a function of the correction to the time derivative 
of the acceleration. 


0 * F(R,V, A,t) 

WHERE F IS AN n ELEMENT VECTOR SUM OF THE STATIC AND DYNAMIC FORCES 
THEN 


0 


- F ^R (0) + 1 aAAt 3 , V (0> + I PAAt 2 , A (0) + AAt, 


OR 


0 = F(AA) 
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LINEARIZED EQUATIONS OF MOTION 

To solve for the correction, the equations of motion are linearized about the 
estimated values. The instantaneous stiffness, damping, and mass are defined 
by the various partial derivatives with respect to displacement, velocity, and 
acceleration. If the linearization is done numerically, the stiffness, damp- 
ing, and mass do not have to be calculated. The numerical differentiation of 
the correction to the time derivative of the acceleration is all that is 
needed. 

This solution procedure is equivalent to the Newton-Raphson technique. The 
numerical differentiation and the solution of the linearized equations of 
motion are computationally time consuming, although straight forward. The 
multigrid technique could be potentially orders of magnitudes faster. The 
linearized equations of motion will be the basis for generating a coarse model 
from a fine model. 


WHERE 


0 - F (0) - BAA 


F (01 = F(R (0> ,V (0) ,A (0) .t) 
B =lot 3 K + Ipt 2 C + tM 


K - 




.-dF 


dA 


BAA - F (0) 
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STRUCTURAL CONDENSATION 

If the linearized equation set is partitioned into nodes belonging to a coarse 
model (upper partition) and the nodes that are eliminated from the fine model 
(lower partition), then structural condensation can be used to solve for the 
coarse model. In addition, the structural condensation process can be used to 
interpolate the solution from the coarse model to the fine model. If the higher 
frequency part of the solution is found on the fine model and the lower fre- 
quency part of the solution is found on the coarse model, then the resultant 
forces must be zero. Thus the solution for the nodes eliminated from the fine 
model can be found. 


l 

• 


(f) 

B 11 ! B 12 

AAj 


_ f l 

\ “■ 
t 

• 


(f) 

B 21 ; b 22 

AA 2 


f 2 J 


( B ll ' B 12 B 22 B 2l) A ^1 = F i f> ‘ B 12 B 22 F 2 f> 

• -1/ (f) . \ 

AA 2 = B 22 \ f 2 ” ^21 AA 1/ 

LET B* c * = ■ ^2^22^21’ F ^ = F 1 ” ^12^22 F 2 

/. B (c >AAj = F (c> 

(f) • -1 

IF F 2 =0 => AA 2 = -B 22 B2iAA 1 (INTERPOLATOR) 
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FINE-TO-COARSE AND COARSE-TO-FINE MODEL TRANSFORMATIONS 


The f ine-to-coarse model transformation is a rectangular matrix that averages 
the force from the fine model to the coarse model. The upper partition is an 
identity matrix, and the lower partition is defined in the structural condensa- 
tion process. The coarse-to-f ine transformation interpolates the correction of 
the time derivative of the acceleration from the coarse to fine model. In the 
symmetric case, the f ine-to-coarse transformation is the transpose of the 
coarse-to-f ine transformation. 


FINE-TO-COARSE MODEL TRANSFORMATION 


$ = 


I 


-1 

’ b 22 b 21 


AA 


COARSE-TO-FINE MODEL TRANSFORMATION 

0 = [i ; " B 12 B 22] =>F (c) =6F (f) 

B AA - F (f) => BBSA/^ = 9F (f) 


OR 


B (c) AA 1 - F (c) 
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NONLINEAR CONDENSATION 


The nonlinear condensation process transforms the independent variables from 
the coarse model to the fine model and the dependent variables from the fine 
to coarse model. Thus the resultant forces are relaxed on the coarse model. 
This would only require the inversion of a diagonal matrix. The corrections 
on the coarse model are then interpolated to the fine model. The lineariza- 
tion of the equations of motion are not needed in the solution process, but 
are needed only to define the transformations. 


0 « F(AA) => 0 - eRSAAj) 


RELAXATION 


0 « 9F (f) - DAAj =* 
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LOCAL STRUCTURAL CONDENSATION 


The linearization of the equations of motion and the structural condensation 
process requires a considerable amount of computational time. Multigrid via 
relaxation is most efficient on nearest neighbor interactions. Thus only a 
partial linearization of the equations of motion is necessary. The equations 
of motion have to be linearized only with respect to the node under considera- 
tion and its nearest neighbors. Applying condensation to this local inter- 
action model results in local structural condensation. In the case of a beam, 
this linearization results in a block tridiagonal matrix and the structural 
condensation results in a coarse model in which every other node is removed 
from the fine model. 




B IS BLOCK TRIDIAGONAL— INCLUDES NEAREST NEIGHBOR INTERACTION, 
NEGLECTS FAR INTERACTION 
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BEAM EXAMPLE 


If the tridiagonal equation set is reordered so that the even numbers are in 
the upper partition for both the fine and coarse model and the odd numbers are 
in the lower partition for the fine model, then the structural condensation has 
a simple form. In the reordered equation set, the block matrices on the diago- 
nal of the partitions are diagonal. The inversions of these block matrices 
are trivial. 
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SOLUTION OF BEAM EXAMPLE 


The solution for the nonidentity partition of both transformations is tridiag- 
onal. For the f ine-to-coarse transformation, the nonidentity partition is also 
lower triangular. For the coarse-to-f ine transformation, the nonidentity par- 
tition is also upper triangular. 



-1 

T L *' Z 2L-1 W 2L-1 


-1 

S L = -Z 2L-1 U 2L-1 



X L - -U 2L Z 2L-1 


Y L = " W 2L Z 2L+1 



.(c) 


Z L =Z 2L + X L W 2L-1 + Y L U 2L+1 
uf -Xl u 2L-1 


(c) 

W L =Y l W 2l+ i 


-(c) 


•(c) 


(c) 

F L = F 2L + X L f 2L-1 + Y L F 2L+1 
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ACCELERATION PARAMETER 


Normally the relaxation technique can be improved by using a weighted average 
of the previous and present calculated values of the corrections to the solu- 
tion (overrelaxation). The rate of convergence of the high-frequency compo- 
nents can be improved at the expense of the low-frequency components. For 
this improvement, an estimate of the highest frequency eigenvalue is needed. 
The Rayleigh quotient is a good method to estimate the highest eigenvalue (at 
least in the symmetric case). In addition the highest eigenvalue should be a 
strong function of the nearest neighbors, therefore local linearization could 
be used in the Rayleigh quotient. 


t(A) IS BASED ON LOCAL COEFFICIENTS 

(Ag) T D(AA~) 

(aA) t (a A) 

THIS IS THE RAYLEIGH QUOTIENT 
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MULTIGRID METHOD 


In summary, the multigrid method for structural dynamics is performed as fol- 
lows. First relax the equations of motion on the fine grid to obtain the high- 
frequency components of the solution. Then calculate the norm of the residual 
on the fine model. Next check to see that the norm is small enough for a 
solution. If not, check to see if the norm has changed significantly from the 
previous iteration. If the norm has changed, then relax the solution until 
the norm stops changing. This indicates that the high-frequency components on 
this model have been found. 

To find the lower frequency components of the solution, use the local struc- 
tural condensation to generate a coarse model. On the coarse model, use relax- 
ation to generate the lower frequency components of the solution. These lower 
frequency components are interpolated to the fine grid where the norm of the 
residual is calculated. Based on this norm, either a solution is found, more 
relaxation is needed, or a coarser model is needed. The process is repeated 
until a solution is found. 


• RELAX ON FINE GRID TO GET HIGH-FREQUENCY COMPONENT 

• CALCULATE RESIDUAL ON FINE GRID 

• CHECK RESIDUAL FOR SOLUTION 

• CHECK CHANGE IN RESIDUAL FOR CHANGE IN GRID 

• STATIC CONDENSE TO COARSE GRID 

• RELAX ON COARSE GRID TO GET LOW-FREQUENCY COMPONENT 

• INTERPOLATE LOW-FREQUENCY TO FINE GRID 
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PHYSICAL MODEL 



FINE MODEL 
COARSE MODEL 

^ f ■ # ■VA f # ' Wr < ♦ 'W e ♦ 

^ VA ♦ ■■■ — ■» ■ m VW 

^ VW - • — 
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MULTIGRID ANALYSIS APPLIED TO 
TRANSMISSION DYNAMICS 


Complete transmission dynamic analyses are rare in the open literature. David 
and Mitchell (1986) have used a modal balance technique. The problem with 
modal techniques is that the nonlinearities cause the set of modes not to be 
closed. This results in side bands around the tooth passing frequency. There- 
fore, the solutions may not always include all of the important modes. Also, 
superfluous modes tend to overwhelm the solution technique. The time march 
multigrid method should eliminate these problems. 
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TRANSMISSION DYNAMIC ANALYSIS 


Transmission dynamics is a case of nonlinear structural dynamics. Physically 
a transmission is composed of gears, shafts, bearings, seals, and a case. The 
case and the shafts can be modeled by finite element methods. The bearings and 
seals are modeled by special programs developed in tribology and other areas. 
Gear interactions are developed for some kinds of gears, but not for others. 
Thus a transmission can be modeled by a number of linear and nonlinear finite 
elements. As a first approximation, a transmission can be modeled as a beam 
structure. The transmission can be analyzed as a multishaft, combined lateral, 
torsional, and axial rotordynamic system. 
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PHYSICAL MODEL 
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SPECIAL FEATURES 


Special features complicate the dynamic analysis of transmissions. Gyroscopic 
and gear forces cause nonlinear lateral and torsional coupling. Gear-tooth 
passing frequencies are high-frequency forcing functions and, therefore, imply 
a need for a fine structural model. Gear-gear interactions cause the system 
to have a wide bandwidth. 


♦ GYROSCOPIC AND GEAR FORCES CAUSE NONLINEAR LATERAL AND TORSIONAL 
COUPLING 


• GEAR-TOOTH PASSING FREQUENCIES ARE HIGH-FREQUENCY FORCING FUNCTIONS- 
IMPLIES NEED FOR FINE STRUCTURAL MODEL 

• GEAR-GEAR INTERACTIONS CAUSE WIDE BAND WIDTH 
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GEAR TOOTH INTERACTION 


Consider gear tooth interaction. For any gear set, the line of force does not 
pass through the gear centers. In the case of spur gears any perturbation of 
the radial force will result in a perturbation of the tangential force and vice 
versa. In the case of helical or spiral gears any perturbation of the radial 
force will result in perturbations of both the axial and tangential force. 

These perturbations result in a nonlinear coupling between the axial, tangen- 
tial, and radial directions. 


SPUR GEAR HELICAL GEAR 
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GEAR TOOTH MODEL 


Consider the gear tooth interactions of a spur gear set* The contact point 
varies as the angle of the gear set varies. Machining errors cause the contact 
point to move. High torque can cause the teeth to bend. The number of teeth 
in contact varies as the torque varies. Negative torque can result in back- 
lash. The force must be transmitted through the contact point. All these 
effects cause nonlinear time varying interactions between the spur gears set. 
For the other kind of gears the interaction is more complicated. Thus, gear 
tooth interactions cause high-frequency forcing functions on the structure. 


r LINE OF ACTION 



CD-88-32955 


I 




1-86 


TYPICAL TRANSMISSION 


In a typical transmission there are many gear sets. Each of these gear sets 
causes one location on the structure to interact with another point on the 
structure. Thus, far interactions are important and the structural model has 
a wide bandwidth. 
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POTENTIAL IMPACT 


The multigrid method, although used for years in fluid dynamics, now offers a 
new approach to nonlinear structural dynamics. The computing time does not 
depend on the cube of the number of degrees of freedom. Thus, dramatic reduc- 
tions in computing time are possible. In addition, the relaxation process is 
applicable to parallel computation. Thus, the method is very attractive for 
future computers. 


• NEW APPROACH TO NONLINEAR STRUCTURAL DYNAMIC SIMULATION 


• DRAMATIC REDUCTION IN COMPUTING TIME 


• APPLICABLE TO PARALLEL COMPUTERS 
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APPENDIX - SYMBOLS 


A acceleration vector 

B linearized coefficient matrix 

C damping matrix 

D diagonal matrix 

F force vector 

G gear torque 

K stiffness matrix 

M mass matrix 

n number of degrees of freedom 

R displacement vector 

r pitch radius of gear 

S block matrix used in 4> (beam solution) 

T block matrix used in 4> (beam solution) 

t time 

U block matrix on lower diagonal of B 

V velocity vector 

W block matrix on upper diagonal of B 
X block matrix used in © (beam solution) 

Y block matrix used in 0 (beam solution) 

Z block matrix on diagonal of B 

a constant, modifying Taylor series 

3 constant, modifying Taylor series 

e weighting factor used in overrelaxation 

0 coarse-to-f ine transformation 

9 helical gear angle 

X highest eigenvalue 

£ value between 0 and t 
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$ f ine-to-coarse transformation 

^ pressure angle 

Subscripts : 

a axial 

L node number 

r radial 

t tangential 

Superscripts : 

c coarse model 

f fine model 

T transpose 

• time derivative 

(0) estimated value 
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